Untitled
Неделя один
Обозначения
Одна зависимая (объясняемая) переменная: y
Несколько регрессоров (предикторов, объясняющих переменных): x, z…
По каждой переменной n наблюдений: \(y_1, y_2...y_n\)
Модель — это формуля для объясняемой переменной.
Пример
Возьмём например данные по машинам 1920-х годов. Тут видимая линейная взаимосвязь.
Данные по машинам 1920-х гдов
Модель для этих данных может иметь вид \(y_i = \beta_1 + \beta_2x_i + \varepsilon_i\)
Что здесь есть что? У нас есть наблюдаемые переменные,
y — это длина тормозного пути
x — это скорость, с которой ехала машина.
Есть неизвестные коэффициенты $ _1, _2$
случайная составляющая, ошибка
То есть, \(\beta_2\) показывает — насколько увеличивается тормозной путь, если машина разгонится на один лишний километр в час. И есть некая случайная составляющая \(ε_i\), это может быть все что угодно:
- водитель по-другому нажимал на тормоз,
- что-то там на дороге было другое,
то есть это та часть, по которой у нас нет возможности предсказать, но, тем не менее, вот эта случайная ошибка она входит в y.
План действий
Придумать адекватную модель
Получить оценки неизвестных параметров \(\widehat\beta_1, \widehat\beta_2\)
Спрогнозировать, заменив неизвестные параметры на оценки \(\widehat y_i = \widehat\beta_1 + \widehat\beta_2 x_i\)
Метод наименьших квадратов
Как найти \(\beta_1, \beta_2\)? Собственно методом наименьших квадратов.
Если мы придумали какие-то оценки, \(\beta_1, \beta_2\) то соответственно возникает такое понятие как ошибка прогноза
\[\widehat \varepsilon_i = y_i - \widehat y_i\]
Есть суммарная ошибка, чтобы суммарная ошибка не занулялась одна в плюс, другая в минус, не компенсировали друг друга, мы возведем в квадрат. И посчитаем сумму квадратов ошибок прогноза, то есть сумма \(\widehat \varepsilon_i ^2\)
\[Q(\widehat\beta_1, \widehat\beta_2) = \sum_{i=1}^{n} \widehat\varepsilon^{2}_{i} = \sum_{i=1}^{n} (y_i - \widehat y_i)^2 \]
Суть метода МНК
Возьмите в качестве оценок такие коэффициенты \(\widehat\beta_1, \widehat\beta_2\) при которых сумма квадратов ошибок прогноза будет минимальна
Случай для одного регрессора
Всё сводится к тому, что в модели \[M_1: y_i = \beta + \varepsilon_i\]
\[\widehat\beta = \bar y\]
Случай для одного регрессора
Случай для двух регрессоров
Тут несколько сложнее, подробности на скриншоте
Случай для двух регрессоров
Случай для множества регрессоров
Что мы получили:
Случай для двух регрессоров
Ещё раз пробежимся по терминам
Случай для двух регрессоров продолжение
Графическое представление
Изобразим на графике, всё то, о чем мы только что проговорили.
Графический вывод
В частности \(\widehat\beta_1\), т.е. точка на оси ординат от пересечении с прямой регрессии, называют пересечением или Intercept-ом
Метод наименьших квадратов подбирает прямую так, чтобы суммарные расстояния от точек до прямой были минимальными.
Графический вывод продолжение
Множественные регрессоры
Случай множественных регрессоров принципиально не отличается от двух регрессоров. Поэтому рассмотрим на примере трёх предикторов.
Множественные регрессоры
Суммы квадратов
Суммы квадратов
Что есть что:
RSS (Residual Sum of Squares) — этот показатель меряет, насколько велики эпсилон с крышкой, насколько они далеко лежат от нуля.
TSS (Total Sum of Squares) — этот показатель меряет, насколько каждый из \(y_i\) не похож на \(y\) среднее. Если \(y_i\) далеко лежит от \(y\) среднее, соответственно, это слагаемое в сумме квадратов будет большим. И вся сумма квадратов будет большой.
ESS (Explained Sum of Squares) — она показывает, насколько прогнозное значение \(y_i\) с крышкой далеко легло от \(y\) среднего.
Лекбез по линейной алгебре
Язык эконометрики во многом — это язык линейной алгебры, нужно его знать.
Обозначения
Обозначения
Маленькой буквой \(y\) мы будем обозначать вектор, то есть столбик из всех игреков, записанных друг под другом: у_1, у_2 и так далее, у_n.
Ну, соответственно, \(х\) маленькое — это все иксы: х_1, х_2 и так далее, х_n, записанные друг под другом.
Аналогично \(\widehat\beta\).
И еще введем обозначение: единичку со стрелочкой. Это будет вектор-столбец, то есть столбик из единичек в количестве n штук, потому что у нас в модели будет n наблюдений.
Тогда для нашей модели
\[\widehat y = \widehat\beta_1 \cdot \overrightarrow 1 + \widehat\beta_2 \cdot x + \widehat\beta_3 \cdot z\]
- Большая буква X — это все регрессоры, помещенные в одну большую табличку чисел, которые называются матрицей.
Таблица регрессоров
- Рассмотрим ещё такое понятие как длина вектора
Длина вектора
Где у нас бывают длины векторов:
Пример длины вектора
Есть одно замечательное применение. Если скалярное произведение векторов равно нулю, значит они перпендикулярны.
Скалярное произведение векторов
x <- c(1,2, -3)
y <- c(-6,0,-2)
sum(x * y)## [1] 0Линейная модель регрессии в n-мерном пространстве
n-мерное простарнство
Есть вектор y, есть вектор из одних 1, я продолжаю этот вектор до прямой и оказывается, что
- в простой модели без регрессоров спроецировав вектор y на эту прямую, я получу вектор y с крышкой — предсказанные значения.
Соответственно, мы получили попутно еще один факт: если любой вектор проецировать на вектор из 1, то получится вектор средних значений
Геометрия множественной регрессии
Напомню, что мы вывели условие первого порядка
Геометрическая интерпретация условий первого порядка
Это означает, что вектора перпендикулярны.
Облачко это все те вектора, которые можно получить, складывая с некоторыми весами вектор x, вектор z и вектор из единичек. Это гиперплоскость.
Мы получаем следующую интрпретацию метода наименьших квадратов:
Первый геометрический факт Прогнозы, которые мы получаем по методу наименьших квадратов — это проекция исходного вектора зависимых переменных y на множество векторов, получаемых с помощью сложения с разными весами вектора из единичек, вектора x и вектора z.
Второй геометрический факт Если я спроецировать вектор y на прямую, порождаемую вектором из единичек, и спроецировать y с крышечкой на ту же самую прямую, то эти проекции попадут в одну и ту же точку.
\(TSS = ESS + RSS\)
\(\frac{ESS}{TSS}= \frac{BC^2}{AB^2}=(\frac{BC}{AB})^2= (cos \varphi)^2\)
Геометрчиеские факты
Коэффициент детерминации
Если в регрессию включён свободный член \(\beta_1\) то действуют следующие правила:
Правила при включении свободного члена в регрессию
Эти правила позволяют придумать простой показатель качества работы модели — коэффициент детерминации.
Чем прогнозы точнее похожи на настоящие значения, тем меньше будут ошибки прогнозов и тем меньше будет сумма квадратов ошибок прогнозов RSS. Соответственно \(\frac{ESS}{TSS} = R^2\) будет примерно равен 1 если RSS будет у ноля. Соответственно мы получим коэффициент детерминации, который всегда лежит от 0 до 1.
Другими словами — коэффициент детерминации это доля объяснённого разброса в общем разбросе.
.
С одной стороны, коэффициент детерминации — это доля объясненной дисперсии игрек, доля объясненного разброса.
С другой стороны, коэффициент детерминации — это выборочная корреляция между прогнозами и настоящим игрек, взятое в квадрат.
Чем коэффициент детерминации выше, тем больше предсказание похожи на реальные значения. Чем коэффициент детерминации выше, тем выше доля объясненной дисперсии.
Пример с фертильностью. МНК в R
Посмотрим как зависит фертильность от других
t <- swiss
# нарисовать много диаграм рассеивания
# install.packages("GGally") # матрица диаграмм рассеяния
library("GGally")## Loading required package: ggplot2str(t)## 'data.frame': 47 obs. of 6 variables:
## $ Fertility : num 80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
## $ Agriculture : num 17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
## $ Examination : int 15 6 5 12 17 9 16 14 12 16 ...
## $ Education : int 12 9 5 7 15 7 7 8 7 13 ...
## $ Catholic : num 9.96 84.84 93.4 33.77 5.16 ...
## $ Infant.Mortality: num 22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...ggpairs(t)
Уже тут можно посмотреть много всяких корреляций.
Перейдём к оценке
model2 <- lm(data = t, Fertility ~ Agriculture + Education + Catholic)
# КОэффициенты
coef(model2)## (Intercept) Agriculture Education Catholic
## 86.2250198 -0.2030377 -1.0721468 0.1452013# Спрогнозированные значения
fitted(model2)## Courtelary Delemont Franches-Mnt Moutier Neuveville
## 71.35382 79.73758 86.36550 76.21257 62.05992
## Porrentruy Broye Glane Gruyere Sarine
## 84.70365 77.94869 77.98965 82.07990 76.37831
## Veveyse Aigle Aubonne Avenches Cossonay
## 81.01451 62.00804 65.34456 61.67811 67.20324
## Echallens Grandson Lausanne La Vallee Lavaux
## 72.85406 71.22373 54.02437 62.00809 62.16632
## Morges Moudon Nyone Orbe Oron
## 64.12130 72.47751 65.22299 69.41765 71.04507
## Payerne Paysd'enhaut Rolle Vevey Yverdon
## 66.61076 70.48740 64.27982 63.09324 68.48321
## Conthey Entremont Herens Martigwy Monthey
## 81.11782 77.02791 80.38838 78.28372 84.09311
## St Maurice Sierre Sion Boudry La Chauxdfnd
## 75.54878 80.27332 73.53528 66.37864 74.87034
## Le Locle Neuchatel Val de Ruz ValdeTravers V. De Geneve
## 70.52554 50.79966 71.80743 76.17918 35.30542
## Rive Droite Rive Gauche
## 52.99371 57.97821# Остатки
residuals(model2)## Courtelary Delemont Franches-Mnt Moutier Neuveville
## 8.8461774 3.3624191 6.1345049 9.5874340 14.8400828
## Porrentruy Broye Glane Gruyere Sarine
## -8.6036475 5.8513084 14.4103471 0.3201012 6.5216935
## Veveyse Aigle Aubonne Avenches Cossonay
## 6.0854871 2.0919628 1.5554442 7.2218873 -5.5032423
## Echallens Grandson Lausanne La Vallee Lavaux
## -4.5540632 0.4762715 1.6756342 -7.7080933 2.9336804
## Morges Moudon Nyone Orbe Oron
## 1.3786987 -7.4775133 -8.6229883 -12.0176461 1.4549265
## Payerne Paysd'enhaut Rolle Vevey Yverdon
## 7.5892409 1.5125978 -3.7798150 -4.7932370 -3.0832083
## Conthey Entremont Herens Martigwy Monthey
## -5.6178154 -7.7279097 -3.0883806 -7.7837172 -4.6931098
## St Maurice Sierre Sion Boudry La Chauxdfnd
## -10.5487835 11.9266828 5.7647206 4.0213575 -9.1703410
## Le Locle Neuchatel Val de Ruz ValdeTravers V. De Geneve
## 2.1744592 13.6003353 5.7925740 -8.5791790 -0.3054172
## Rive Droite Rive Gauche
## -8.2937095 -15.1782122# ПОказатель RSS
deviance(model2)## [1] 2567.884# Показатель R квадрат
report <- summary(model2)
report$r.squared## [1] 0.6422541# Это тоже самое что
cor(t$Fertility, fitted(model2))^2## [1] 0.6422541Тесты
data("mtcars")
m1 <- lm(data = mtcars, mpg ~ disp + hp + wt)
summary(m1)$r.squared## [1] 0.8268361m2 <- lm(data = mtcars, mpg ~ cyl + hp + wt)
summary(m2)$r.squared## [1] 0.84315m3 <- lm(data = mtcars, mpg ~ disp + cyl + wt)
summary(m3)$r.squared## [1] 0.832607m4 <- lm(data = mtcars, mpg ~ disp + hp + cyl)
summary(m4)$r.squared## [1] 0.7678877m5 <- lm(data = mtcars, mpg ~ disp + hp + wt + am)
coef(m5)## (Intercept) disp hp wt am
## 34.209443370 0.002489354 -0.039323213 -3.046747000 2.159270737m6 <- lm(data = mtcars, mpg ~ disp + hp + wt)
summary(m6)$r.squared## [1] 0.8268361data("sleep")
var(tail(sleep, 11)$extra)## [1] 3.618max(sleep$extra) - min(sleep$extra)## [1] 7.1mean(sleep$extra)^3## [1] 3.652264sleep[8,]## extra group ID
## 8 0.8 1 8